#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "design.h"

void setup_simplified (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k;

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 0)
			continue;
        
		simplified[j][2 * N + 1] = simplified[j][2 * N + 2] = 0;
		for (k = n; k < N; k++)
		{
			simplified[j][2 * N + 1] += simplified[j][k];
			simplified[j][2 * N + 2] += simplified[j][N + k];
		}
	}
}

void apply_Fr (int N, int L, int simplified[L][2 * N + 3], int Fr[N], int n)
{
	int j, k;

	/* Apply Forbidden R */
	for (k = 0; k < N; k++)
	{
		if (Fr[k] == 1)
		{
			for (j = 0; j < L; j++)
			{
				simplified[j][N + k] = 0;
			}
		}
	}

	/* Setup Array simplified */
	setup_simplified (N, L, simplified, n);
}

void apply_Fr2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int j;

	/* Apply Forbidden R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == -1 && simplified[j][N + p] == 1)
		{
			simplified[j][N + p] = 0;
			simplified[j][2 * N + 2] --;
		}
	}
}

void simplify_Fr (int N, int L, int state[L][2 * N + 1], int simplified[L][2 * N + 3], int Fr[N], int i, int mode, int mode2)
{
	int j, k, flag;

	/* Get Forbidden R */
	for (j = 0; j < L; j++)
	{
		flag = 0;

		if (j == 0 && mode == 2)
		{
			flag = 1;
		}
		else if (mode2 == 0)
		{
			flag = 1;
		}
		else if (mode2 == 1)
		{
			if (state[j][i] == state[j][N + i])
				flag = 1;
		}
		else if (mode2 == 2)
		{
			if (state[j][N + i] != state[j][2 * N])
				flag = 1;
		}
		else if (mode2 == 3)
		{
			if (state[j][i] == state[j][N + i] && state[j][N + i] != state[j][2 * N])
				flag = 1;
		}

		if (state[j][N + i] == 0)
		{
			simplified[j][2 * N] = -1;
			if (flag)
			{
				for (k = 0; k < N; k++)
				{
					if (state[j][k] == 1)
					{
						simplified[j][k] = 1;
						if (i != k)
							simplified[j][N + k] = 1;
					}
				}
			}
			else
			{
				for (k = 0; k < N; k++)
				{
					if (i == k)
					{
						if (state[j][k] == 1)
							simplified[j][k] = 1;
					}
					else
					{
						if (state[j][k] == 1 || state[j][N + k] == 1)
							simplified[j][k] = 1;
                                                
						if (state[j][k] == 1 && state[j][N + k] == 1)
							simplified[j][N + k] = 1;
					}
				}
			}
		}
		else
		{
			simplified[j][2 * N] = 1;
			if (flag)
			{
				for (k = 0; k < N; k++)
				{
					if (state[j][k] == 1)
					{
						simplified[j][k] = 1;
						if (i != k)
							Fr[k] = 1;
					}
				}
			}
			else
			{
				for (k = 0; k < N; k++)
				{
					if (i == k)
					{
						if (state[j][k] == 1)
							simplified[j][k] = 1;
					}
					else
					{
						if (state[j][k] == 1 && state[j][N + k] == 1)
							simplified[j][k] = 1;

						if (state[j][k] == 1 || state[j][N + k] == 1)
							Fr[k] = 1;
					}
				}
			}
		}
	}

	/* Apply Forbidden R and Reset Array simplified */
	apply_Fr (N, L, simplified, Fr, 0);
}

void apply_Fg (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int Fg[N], int i, int n)
{
	int j, k;

	/* Apply Forbidden G */
	for (k = n; k < N; k++)
	{
		if ((i == k && Rr[k] == 1) || (i != k && Fg[k] == 1))
		{
			for (j = 0; j < L; j++)
			{
				if (simplified[j][k] == 1)
				{
					simplified[j][k] = 0;
					simplified[j][2 * N + 1] --;
					if (simplified[j][2 * N] == -1 && simplified[j][2 * N + 1] == 0)
					{
						simplified[j][2 * N] = 0;
					}
				}
			}
		}
	}
}
		
void apply_Fg2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int k;

	/* Apply Forbidden G */
	for (k = 0; k < L; k++)
	{
		if (simplified[k][2 * N] != 0 && simplified[k][p] == 1)
		{
			simplified[k][p] = 0;
			simplified[k][2 * N + 1] --;
			if (simplified[k][2 * N] == -1 && simplified[k][2 * N + 1] == 0)
				simplified[k][2 * N] = 0;
		}
	}
}

/* Return 1: find new member in Forbidden G */
int simplify_Fg (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int Fg[N], int i, int n)
{
	int j, k, flag;

	flag = 0;

	/* Get Forbidden G */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1 || simplified[j][2 * N + 2] != 0)
			continue;

		flag = 1;
		simplified[j][2 * N] = 0;
		for (k = n; k < N; k++)
		{
			if (simplified[j][k] == 1)
			{
				if (i == k)
					Rr[k] = 1;
				else
					Fg[k] = 1;
			}
		}
	}

	return flag;
}

void apply_Rg (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Fr[N], int i, int n)
{
	int j, k, p;

	/* Apply Required G */
	for (k = n; k < N; k++)
	{
		if ((i == k && Fr[k] == 0) || (i != k && Rg[k] == 0))
			continue;

		for (j = 0; j < L; j++)
		{
			if (simplified[j][k] == 1)
			{
				if (simplified[j][2 * N] == 1)
				{
					simplified[j][2 * N] = 0;
				}
				else if (simplified[j][2 * N] == -1)
				{
					for (p = n; p < N; p++)
					{
						simplified[j][p] = 0;
					}
					simplified[j][2 * N + 1] = 0;
				}
			}
		}
	}
}

void apply_Rg2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int j, k;

	/* Apply Required R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != 0 && simplified[j][p] == 1)
		{
			if (simplified[j][2 * N] == 1)
			{
				simplified[j][2 * N] = 0;
			}
			else
			{
				simplified[j][2 * N + 1] = 0;
				for (k = 0; k < N; k++)
				{
					simplified[j][k] = 0;
				}
			}
		}
	}
}

/* Return 1: find new member in Required G */
int simplify_Rg (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Fr[N], int i, int n)
{
	int j, k, p, flag;

	flag = 0;

	/* Get Required G */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 1 && simplified[j][2 * N + 1] == 1)
		{
			flag = 1;
			simplified[j][2 * N] = 0;

			for (k = n; k < N; k++)
			{
				if (simplified[j][k] == 1)
				{
					if (i == k)
						Fr[k] = 1;
					else
						Rg[k] = 1;

					break;
				}
			}
		}
	}

	return flag;
}

void combinationR (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k, p, flag;

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1)
			continue;

		for (k = j + 1; k < L; k++)
		{
			if (simplified[k][2 * N] != -1)
				continue;

			if (simplified[j][2 * N + 2] != simplified[k][2 * N + 2])
				continue;

			flag = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[j][N + p] != simplified[k][N + p])
				{
					flag = 1;
					break;
				}
			}

			if (flag == 1)
				continue;

			/* Starting combination */
			simplified[k][2 * N] = 0;

			if (simplified[j][2 * N + 1] == 0)
				continue;

			if (simplified[k][2 * N + 1] == 0)
			{
				simplified[j][2 * N + 1] = 0;

				for (p = n; p < N; p++)
					simplified[j][p] = 0;

				continue;
			}

			simplified[j][2 * N + 1] = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[k][p] == 1)
					simplified[j][p] = 1;
				simplified[j][2 * N + 1] += simplified[j][p];
			}
		}
	}
}

void elimate1 (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k, p, flag;

	/* elimate R equation from G equation */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != 1)
			continue;

		for (k = 0; k < L; k++)
		{
			if (j == k)
				continue;

			if (simplified[k][2 * N] == 0)
				continue;

			if (simplified[k][2 * N + 1] < simplified[j][2 * N + 1])
				continue;

			flag = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[j][p] == 1 && simplified[k][p] == 0)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 0)
			{
				if (simplified[k][2 * N] == 1)
				{
					simplified[k][2 * N] = 0;
				}
				else
				{
					simplified[k][2 * N + 1] = 0;
					for (p = n; p < N; p++)
					{
						simplified[k][p] = 0;
					}
				}
			}
		}
	}
}

void apply_Rr (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int i, int n)
{
	int j, k;

	/* Apply Required R */
	for (k = n; k < N; k++)
	{
		if (i == k || Rr[k] == 0)
			continue;

		for (j = 0; j < L; j++)
		{
			if (simplified[j][2 * N] == -1 && simplified[j][N + k] == 1)
			{
				simplified[j][2 * N] = 0;
			}
		}
	}
}

void apply_Rr2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int j;

	/* Apply Required R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == -1 && simplified[j][N + p] == 1)
		{
			simplified[j][2 * N] = 0;
		}
	}
}

/* Return 1: find new member in Required G ; Return 2: find new member in Forbiddgen R ; Return 3: find both */
int simplify_Rr (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int Fg[N], int n)
{
	int j, k, p, flag1, flag2;

	flag1 = flag2 = 0;

	/* Find Required R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1)
			continue;

		if (simplified[j][2 * N + 2] == 1)
		{
			for (k = n; k < N; k++)
			{
				if (simplified[j][N + k] == 1)
					break;
			}

			p = k;

			if (simplified[j][2 * N + 1] == 0)
			{
				flag1 = 1;
				Rr[p] = 1;
				simplified[j][2 * N] = 0;
			}

			if (Fg[p] == 0)
			{
				flag2 = 1;
				Fg[p] = 1;

				/* Apply Forbidden G */
				apply_Fg2 (N, L, simplified, p);
			}
		}
	}

	if (flag1 == 1 && flag2 == 0)
		return 1;
	else if (flag1 == 0 && flag2 == 1)
		return 2;
	else if (flag1 == 1 && flag2 == 1)
		return 3;
	else
		return 0;
}

void elimate2 (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k, p, flag;

	/* elimate R equation from R equation */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1)
			continue;

		for (k = 0; k < L; k++)
		{
			if (j == k)
				continue;

			if (simplified[k][2 * N] != -1)
				continue;

			if (simplified[j][2 * N + 2] > simplified[k][2 * N + 2])
				continue;

			if (simplified[j][2 * N + 1] == 0)
				;
			else if (simplified[k][2 * N + 1] == 0)
				continue;
			else if (simplified[j][2 * N + 1] < simplified[k][2 * N + 1])
				continue;

			flag = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[j][2 * N + 1] != 0 && simplified[k][p] == 1 && simplified[j][p] == 0)
				{
					flag = 1;
					break;
				}

				if (simplified[k][N + p] == 0 && simplified[j][N + p] == 1)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 0)
			{
				simplified[k][2 * N] = 0;
			}
		}
	}
}

int get_optional (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int i, int n)
{
	int j, k, flag;

	for (j = n; j < N; j++)
	{
		Og[j] = Or[j] = 0;
	}

	/* Get the option G */
	for (j = n; j < N; j++)
	{
		if (Rg[j] == 1 || Fg[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] == 0)
				continue;

			if (simplified[k][j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
			Og[j] = 1;
	}

	/* Get the option R */
	for (j = n; j < N; j++)
	{
		if (i == j || Rr[j] == 1 || Fr[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] != -1)
				continue;

			if (simplified[k][N + j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
			Or[j] = 1;
	}

	/* self inhibition */
	if (Rr[i] == 0 && Fr[i] == 0 && n <= i)
	{
		flag = 0;
		for (j = 0; j < L; j++)
		{
			if (simplified[j][2 * N] == 0)
				continue;

			if (simplified[j][i] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
			Or[i] = 1;
	}

	return 0;
}

int get_removable (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int i, int n)
{
	int j, k, flag, p;

	p = 0;

	/* Get the removable G */
	for (j = n; j < N; j++)
	{
		if (Rg[j] == 1 || Fg[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] == 1 && simplified[k][j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
		{
			p = Fg[j] = 1;
			apply_Fg2 (N, L, simplified, j);
		}
	}

	/* Get the removable R */
	for (j = n; j < N; j++)
	{
		if (i == j || Rr[j] == 1 || Fr[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] == -1 && simplified[k][N + j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
		{
			p = Fr[j] = 1;
			apply_Fr2 (N, L, simplified, j);
		}
	}

	/* self inhibition */
	if (Rr[i] == 0 && Fr[i] == 0 && n <= i)
	{
		flag = 0;
		for (j = 0; j < L; j++)
		{
			if (simplified[j][2 * N] == -1 && simplified[j][i] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
		{
			p = Fr[i] = 1;
			apply_Rg2 (N, L, simplified, i);
		}
	}

	return p;
}

void swap (int *a, int *b)
{
	int tmp;

	tmp = *a;
	*a = *b;
	*b = tmp;
}

int reorder (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int edge[N], int order[N], int i, int j, int p)
{
	int k;

	if (i == j)
		return p;
	else
	{
		for (k = 0; k < L; k++)
		{
			swap (&simplified[k][i], &simplified[k][j]);
			swap (&simplified[k][N + i], &simplified[k][N + j]);
		}

		swap (&Rg[i], &Rg[j]);
		swap (&Rr[i], &Rr[j]);
		swap (&Fg[i], &Fg[j]);
		swap (&Fr[i], &Fr[j]);
		swap (&Og[i], &Og[j]);
		swap (&Or[i], &Or[j]);
		swap (&edge[i], &edge[j]);
		swap (&order[i], &order[j]);
	}

	if (i == p)
		return j;
	else if (j == p)
		return i;
	else
		return p;
}

int check_irreducible (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int edge[N], int i)
{
	int j, k, flag, flag2, p;
	int Rg2[N], Rr2[N];

	for (j = 0; j < N; j++)
	{
		Rg2[j] = Rr2[j] = 0;
	}

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 0)
			break;
		else if (simplified[j][2 * N] == 1)
		{
			if (simplified[j][i] == 1 && edge[i] == 0)
				continue;

			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (k == i)
					continue;
				else if (simplified[j][k] == 1 && edge[k] == 1)
				{
					flag ++;
					if (flag == 1)
						p = k;
					else
						break;
				}
			}

			if (flag == 1)
				Rg2[p] = 1;
		}
		else if (simplified[j][2 * N] == -1)
		{
			if (simplified[j][2 * N + 1] > 0)
			{
				flag = 0;
				for (k = 0; k < N; k++)
				{
					if (k == i)
						continue;
					else if (simplified[j][k] == 1 && edge[k] == 1)
					{
						flag = 1;
						break;
					}
				}
                                
				if (flag == 0)
				{
					if (simplified[j][i] == 0)
						continue;
					else if (edge[i] == -1)
					{
						flag = 1;
						p = i;
					}
				}
				else
					flag = 0;
			}
			else
				flag = 0;

			for (k = 0; k < N; k++)
			{
				if (simplified[j][N + k] == 1 && edge[k] == -1)
				{
					flag ++;
					if (flag == 1)
						p = k;
					else
						break;
				}
			}

			if (flag == 1)
				Rr2[p] = 1;
		}
	}

	flag = 1;
	for (j = 0; j < N; j++)
	{
		if (edge[j] != 0 && Rg[j] + Rr[j] + Rg2[j] + Rr2[j] == 0)
		{
			flag = 0;
			break;
		}
	}

	return flag;
}

/*
void print (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int n, int edge[N], int order[N])
{
	int j, k, order2[N];

	for (j = 0; j < N; j++)
	{
		order2[order[j]] = j;
	}

	for (j = 0; j < N; j++)
	{
		if (order2[j] >= n)
			printf ("- ");
		else if (edge[order2[j]] == 0)
			printf ("n ");
		else if (edge[order2[j]] == 1)
			printf ("g ");
		else
			printf ("r ");
	}
	printf ("\n");

	for (j = 0; j < L; j++)
	{
		for (k = 0; k < 2 * N + 3; k++)
		{
			printf ("%d ", simplified[order2[j]][k]);
		}
		printf ("\n");
	}
	printf ("\n");

	printf ("Rg : ");
	for (j = 0; j < N; j++)
	{
		printf (" %d", Rg[order2[j]]);
	}
	printf ("\n");

	printf ("Rr : ");
	for (j = 0; j < N; j++)
	{
		printf (" %d", Rr[order2[j]]);
	}
	printf ("\n");

	printf ("Fg : ");
	for (j = 0; j < N; j++)
	{
		printf (" %d", Fg[order2[j]]);
	}
	printf ("\n");

	printf ("Fr : ");
	for (j = 0; j < N; j++)
	{
		printf (" %d", Fr[order2[j]]);
	}
	printf ("\n");

	printf ("Og : ");
	for (j = 0; j < N; j++)
	{
		printf (" %d", Og[order2[j]]);
	}
	printf ("\n");

	printf ("Or : ");
	for (j = 0; j < N; j++)
	{
		printf (" %d", Or[order2[j]]);
	}
	printf ("\n");

	printf ("\n");
}
*/

int check_solution (int N, int L, int simplified[L][2 * N + 3], int edge[N], int i, int old, int new)
{
	int j, k, ok, flag, flagR, flagG;

	flagR = flagG = 0;
	if (old == 1)
	{
		flagG = 1;
	}
	else if (old == -1)
	{
		flagR = 1;
	}
	else
	{
		if (new == 1)
			flagR = 1;
		else
			flagG = 1;
	}

	ok = 1;
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 0)
			break;
		else if (flagG && simplified[j][2 * N] == 1)
		{
			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (simplified[j][k] == 1 && k != i && edge[k] == 1)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 0 && (simplified[j][i] == 0 || edge[i] == -1))
			{
				ok = 0;
				break;
			}
		}
		else if (flagR && simplified[j][2 * N] == -1)
		{
			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (simplified[j][N + k] == 1 && edge[k] == -1)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 1)
				continue;

			if (simplified[j][2 * N + 1] == 0 || (simplified[j][i] == 1 && edge[i] == 0))
				flag = 1;
			else
			{
				for (k = 0; k < N; k++)
				{
					if (k == i)
						continue;
                                
					if (simplified[j][k] == 1 && edge[k] == 1)
					{
						flag = 1;
						break;
					}
				}
			}

			if (flag == 1)
			{
				ok = 0;
				break;
			}
		}
	}

	return ok;
}

int check_degree (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int edge[N], int i, int dgflag)
{
	int j, degree, tmp;
	int flag[3]; /* 0-R, 1-N, 2-G */

	degree = 0;

	for (j = 0; j < N; j++)
	{
		if (Rg[j] == 1 || Rr[j] == 1)
			continue;

		flag[0] = flag[1] = flag[2] = -1;

		flag[edge[j] + 1] = 1;

		if (Or[j] == 1)
			flag[0] = flag[1] = 1;
		else if (Fr[j] == 1)
			flag[0] = 0;

		if (Og[j] == 1)
			flag[1] = flag[2] = 1;
		else if (Fg[j] == 1)
			flag[2] = 0;

		if (dgflag != 2)
		{
			if (edge[j] == -1)
				flag[2] = 0;
			else if (edge[j] == 1)
				flag[0] = 0;
		}

		tmp = edge[j];

		if (flag[1] == -1)
		{
			edge[j] = 0;
			flag[1] = check_solution (N, L, simplified, edge, i, tmp, edge[j]);
		}

		if (flag[2] == -1)
		{
			if (flag[1] == 1)
			{
				edge[j] = 1;
				flag[2] = check_solution (N, L, simplified, edge, i, tmp, edge[j]);
			}
			else
			{
				flag[2] = 0;
			}
		}

		if (flag[0] == -1)
		{
			if (j == i && flag[1] == 1)
			{
				edge[j] = -1;
				flag[0] = check_solution (N, L, simplified, edge, i, tmp, edge[j]);
			}
			else
			{
				flag[0] = flag[1];
			}
		}

		edge[j] = tmp;

		flag[edge[j] + 1] = 0;
		degree += flag[0] + flag[1] + flag[2];
	}

	return degree;
}

void init_solution (struct solution *solution)
{
	solution->max = solution->curr = 0;
	solution->edge = NULL;
	solution->disqualify = NULL;
	solution->count = NULL;
}

void free_solution (struct solution *solution)
{
	solution->max = solution->curr = 0;
	free (solution->edge);
}

void get_solution (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int edge[N], int i, int dgflag, int svflag, int pnflag, struct solution *solution, int *degree_dist)
{
	int j;

	if (svflag)
	{
		if (solution->max == solution->curr)
		{
			solution->max += 100;
			solution->edge = realloc (solution->edge, sizeof (int) * N * solution->max);
		}
        
		memcpy (solution->edge + solution->curr * N, edge, sizeof (int) * N);
		solution->curr ++;
	}

	if (pnflag)
	{
		printf ("Node %2d :", i);
		for (j = 0; j < N; j++)
		{
			if (edge[j] == 0)
				printf (" n");
			else if (edge[j] == 1)
				printf (" g");
			else if (edge[j] == -1)
				printf (" r");
		}
		printf ("\n");
	}

	if (dgflag)
	{
		int degree;

		degree = check_degree (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, edge, i, dgflag);
		degree_dist[degree] ++;
	}
}

void get_solution2 (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int edge[N], int i, int dgflag, int svflag, int pnflag, struct solution *solution, int *degree_dist, int n)
{
	if (n == N)
	{
		get_solution (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, edge, i, dgflag, svflag, pnflag, solution, degree_dist);
	}
	else
	{
		if (Og[n] == 1)
		{
			edge[n] = 1;
			get_solution2 (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, edge, i, dgflag, svflag, pnflag, solution, degree_dist, n + 1);
		}

		if (Or[n] == 1)
		{
			edge[n] = -1;
			get_solution2 (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, edge, i, dgflag, svflag, pnflag, solution, degree_dist, n + 1);
		}

		if (Og[n] == 1 || Or[n] == 1)
		{
			edge[n] = 0;
			get_solution2 (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, edge, i, dgflag, svflag, pnflag, solution, degree_dist, n + 1);
		}
		else
		{
			get_solution2 (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, edge, i, dgflag, svflag, pnflag, solution, degree_dist, n + 1);
		}
	}
}

int simplify (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Rr[N], int Fg[N], int Fr[N], int Og[N], int Or[N], int simplified0[L][2 * N + 3], int Rg0[N], int Rr0[N], int i, int n, int edge[N], int flag, int mode, int irflag, int bdflag, int dgflag, int svflag, int pnflag, struct designability *designability, int Link_dist[N + 1 - mode], int degree_dist[2 * (N - mode) + 1], struct solution *solution, int order[N])
{
	int j, k, p;
	int flag1, flag2, flag3;

	int Rg2[N], Rr2[N], Fg2[N], Fr2[N], order2[N];
	int simplified2[L][2 * N + 3];

	memcpy (Rg2, Rg, sizeof (int) * N);
	memcpy (Rr2, Rr, sizeof (int) * N);
	memcpy (Fg2, Fg, sizeof (int) * N);
	memcpy (Fr2, Fr, sizeof (int) * N);
	memcpy (order2, order, sizeof (int) * N);
	memcpy (simplified2, simplified, sizeof (int) * L * (2 * N + 3));

	if (flag == 1)
	{
		edge[n - 1] = 1;
		apply_Rg2 (N, L, simplified2, n - 1);
		if (Fr2[n - 1] == 0)
		{
			apply_Fr2 (N, L, simplified2, n - 1);
		}
	}
	else if (flag == 0)
	{
		edge[n - 1] = 0;
		if (i == n - 1)
		{
			apply_Rg2 (N, L, simplified2, n - 1);
		}
		else
		{
			if (Fr2[n - 1] == 0)
			{
				apply_Fr2 (N, L, simplified2, n - 1);
			}
			if (Fg2[n - 1] == 0)
			{
				apply_Fg2 (N, L, simplified2, n - 1);
			}
		}
	}
	else if (flag == -1)
	{
		edge[n - 1] = -1;
		if (i == n - 1)
		{
			apply_Fg2 (N, L, simplified2, n - 1);
		}
		else
		{
			apply_Rr2 (N, L, simplified2, n - 1);
			if (Fg2[n - 1] == 0)
			{
				apply_Fg2 (N, L, simplified2, n - 1);
			}
		}
	}

	if (mode == 1 && n == 1 && edge[0] == 1 && Fr2[i] == 0)
	{
		Fr2[i] = 1;
		apply_Fr2 (N, L, simplified2, i);
	}

	if (flag != 2)
	{
		flag3 = 0;
                
		do
		{
			p = 0;
			if (flag != -1)
			{
				p = simplify_Fg (N, L, simplified2, Rr2, Fg2, i, n);
				if (p == 1)
					apply_Fg (N, L, simplified2, Rr2, Fg2, i, n);
			}
                        
			flag1 = 0;
			flag2 = 2;
                
			while (flag2 >= 2)
			{
				if (p == 1)
				{
					for (j = 0; j < L; j++)
					{
						if (simplified2[j][2 * N] == 1 && simplified2[j][2 * N + 1] == 0)
							return 0;
					}
				}
                
				if (flag != 1)
				{
					p = simplify_Rg (N, L, simplified2, Rg2, Fr2, i, n);
					if (p == 1)
						apply_Rg (N, L, simplified2, Rg2, Fr2, i, n);
				}
                
				if (flag != -1 || i == n - 1 || flag3 == 1)
				{
					/* Get Required R */
					if (flag1 == 0) /* Combine the R equation */
					{
						combinationR (N, L, simplified2, n);
					}
                                        
					elimate1 (N, L, simplified2, n);
                                        
					flag2 = simplify_Rr (N, L, simplified2, Rr2, Fg2, n);
					if (flag2 == 1 || flag2 ==3)
						apply_Rr (N, L, simplified2, Rr2, i, n);
                                        
                                
					elimate2 (N, L, simplified2, n);
				}
				else
					break;
			}
                
			if (irflag == 0)
				flag1 = get_optional (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, i, n);
			else
				flag1 = get_removable (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, i, n);
			flag3 = 1;
		}
		while (flag1);
	}

	if (n == 0)
	{
		memcpy (Rg0, Rg2, sizeof (int) * N);
		memcpy (Rr0, Rr2, sizeof (int) * N);
		for (j = k = 0; j < L; j++)
		{
			if (simplified2[j][2 * N] != 0)
			{
				memcpy (simplified0[k], simplified2[j], sizeof (int) * (2 * N + 3));
				k ++;
			}
		}
		for (j = k; j < L; j++)
			simplified0[j][2 * N] = 0;
	}

	//print (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, n, edge, order);

	if (n == N)
	{
		int edge2[N];

		/* recover the order */
		for (j = 0; j < N; j++)
		{
			edge2[order[j]] = edge[j];
		}
        
		i = order[i];

		if (irflag && check_irreducible (N, L, simplified0, Rg0, Rr0, edge2, i) == 0)
			return 0;

		p = 0;

		for (j = 0; j < N; j++)
		{
			p += edge[j] * edge[j];
		}

		if (p < designability->min_link[i - mode])
		{
			designability->min_link[i - mode] = p;
		}

		if (bdflag && irflag)
			Link_dist[p] ++;
		else if (bdflag && !irflag)
		{
			int link_dist[N + 1], link_dist1[N + 1], link_dist2[N + 1], q, r;
                        
			link_dist[p] = 1;
			for (j = p + 1; j <= N; j++)
				link_dist[j] = 0;
                        
			for (j = 0; j < N; j++)
			{
				r = Og[j] + Or[j];
				if (r == 1)
				{
					memcpy (link_dist1, link_dist, sizeof (int) * (N + 1));
					link_dist2[p] = 0;
					for (k = p; k < N; k++)
					{
						link_dist2[k + 1] = link_dist1[k];
					}
					for (k = p; k <= N; k++)
					{
						link_dist[k] = link_dist1[k] + link_dist2[k];
					}
				}
				else if (r == 2)
				{
					memcpy (link_dist1, link_dist, sizeof (int) * (N + 1));
					link_dist2[p] = 0;
					for (k = p; k < N; k++)
					{
						link_dist2[k + 1] = 2 * link_dist1[k];
					}
					for (k = p; k <= N; k++)
					{
						link_dist[k] = link_dist1[k] + link_dist2[k];
					}
				}
			}
                        
			for (j = p; j <= N; j++)
				Link_dist[j] += link_dist[j];
		}

		if (dgflag || svflag || pnflag)
		{
			if (irflag)
			{
				get_solution (N, L, NULL, NULL, NULL, NULL, NULL, NULL, NULL, edge2, i, dgflag, svflag, pnflag, solution, NULL);
			}
			else
			{
				int Og3[N], Or3[N];
                                
				for (j = 0; j < N; j++)
				{
					Og3[order[j]] = Og[j];
					Or3[order[j]] = Or[j];
				}

				if (dgflag)
				{
					int Rg3[N], Rr3[N], Fg3[N], Fr3[N];

					for (j = 0; j < N; j++)
					{
						Rg3[order[j]] = Rg2[j];
						Rr3[order[j]] = Rr2[j];
						Fg3[order[j]] = Fg2[j];
						Fr3[order[j]] = Fr2[j];
					}

					get_solution2 (N, L, simplified0, Rg3, Rr3, Fg3, Fr3, Og3, Or3, edge2, i, dgflag, svflag, pnflag, solution, degree_dist, 0);
				}
				else
				{
					get_solution2 (N, L, NULL, NULL, NULL, NULL, NULL, Og3, Or3, edge2, i, dgflag, svflag, pnflag, solution, NULL, 0);
				}
			}
		}

		return 1;
	}

	p = 1;

	if (mode == 1 && n == 0)
	{
		if (Rg2[n] == 1)
		{
			edge[n] = 1;
			return simplify (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, simplified0, Rg0, Rr0, i, n + 1, edge, 2, mode, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution, order2);
		}
	}
	else
	{
		for (j = n; j < N; j++)
		{
			if (Rr2[j] == 1)
			{
				edge[j] = -1;
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, j, i);
				n ++;
			}
			else if (Rg2[j] == 1)
			{
				edge[j] = 1;
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, j, i);
				n ++;
			}
			else if (Fg2[j] == 1 && Fr2[j] == 1)
			{
				edge[j] = 0;
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, j, i);
				n ++;
			}
			else if (irflag)
				;
			else if (Og[j] == 1 && Or[j] == 1)
			{
				edge[j] = 0;
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, j, i);
				n ++;
				p *= 3;
			}
			else if (Og[j] == 1 && Fr2[j] == 1)
			{
				edge[j] = 0;
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, j, i);
				n ++;
				p *= 2;
			}
			else if (Or[j] == 1 && Fg2[j] == 1)
			{
				edge[j] = 0;
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, j, i);
				n ++;
				p *= 2;
			}
		}
        
		if (n == N)
		{
			return p * simplify (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, simplified0, Rg0, Rr0, i, n, edge, 2, mode, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution, order2);
		}
		else if (n == N - 1)
			;
		else
		{
			int countG[N], countR[N], max, max_p;
                        
			for (j = n; j < N; j++)
			{
				countG[j] = countR[j] = 0;
				for (k = 0; k < L; k++)
				{
					if (simplified2[k][2 * N] == 0)
						continue;
					else if (simplified2[k][2 * N] == 1 && simplified2[k][j] == 1)
						countG[j] ++;
					else if (simplified2[k][2 * N] == -1 && simplified2[k][N + j] == 1)
						countR[j] ++;
                        
				}
			}
                        
			max = countR[n];
			max_p = n;
			for (j = n + 1; j < N; j++)
			{
				if (countR[j] > max)
				{
					max = countR[j];
					max_p = j;
				}
			}
			if (max > 0)
			{
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, max_p, i);
			}
			else
			{
				max = countG[n];
				max_p = n;
				for (j = n + 1; j < N; j++)
				{
					if (countR[j] > max)
					{
						max = countG[j];
						max_p = j;
					}
				}
				i = reorder (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, edge, order2, n, max_p, i);
			}
		}
	}

	int NS;

	NS = 0;

	Og[n] = Or[n] = 0;

	if (Fg2[n] == 0)
	{
		NS += simplify (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, simplified0, Rg0, Rr0, i, n + 1, edge, 1, mode, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution, order2);
	}

	NS += simplify (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, simplified0, Rg0, Rr0, i, n + 1, edge, 0, mode, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution, order2);

	if (Fr2[n] == 0)
	{
		NS += simplify (N, L, simplified2, Rg2, Rr2, Fg2, Fr2, Og, Or, simplified0, Rg0, Rr0, i, n + 1, edge, -1, mode, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution, order2);
	}

	return p * NS;
}

int check_satisfy1 (int N, int L, int state[L][2 * N + 1], int i, int mode, int mode1, int mode2, int irflag, int bdflag, int dgflag, int svflag, int pnflag, struct designability *designability, int Link_dist[N + 1 - mode1], int degree_dist[2 * (N - mode1) + 1], struct solution *solution)
{
	int Rg[N], Rr[N]; /* Rg/Rr: Required G/R */
	int Fg[N], Fr[N]; /* Fg/Fr: Forbidden G/R */
	int Og[N], Or[N]; /* Og/Or: Optional G/R */
	int simplified[L][2 * N + 3]; /* [0, N):g; [N, 2*N):r; 2*N:-1/0/1; 2*N+1/2:sum of g/r */
	int simplified0[L][2 * N + 3], Rg0[N], Rr0[N];
	int j, k, p, flag, flag1, flag2;

	/* Set the initial value */
	for (j = 0; j < N; j++)
	{
		Rg[j] = Rr[j] = Fg[j] = Fr[j] = Og[j] = Or[j] = 0;
	}
	Fg[i] = 1; /* No self-activation */
	if (mode1 == 1) Fr[0] = 1;

	if (mode == 2 && state[0][i] == state[0][N + i])        
	{
		Fr[mode1] = Fg[mode1] = 1;         
	}

	for (j = 0; j < L; j++)
	{
		for (k = 0; k < 2 * N + 3; k++)
		{
			simplified[j][k] = 0;
		}
	}

	simplify_Fr (N, L, state, simplified, Fr, i, mode, mode2);

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 1 && simplified[j][2 * N + 1] == 0)
			return 0;
	}

	int edge[N], order[N];

	for (j = 0; j < N; j++)
	{
		order[j] = j;
	}

	designability->NS2[i - mode1] = simplify (N, L, simplified, Rg, Rr, Fg, Fr, Og, Or, simplified0, Rg0, Rr0, i, 0, edge, -2, mode1, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution, order);

	if (designability->NS2[i - mode1] == 0)
		return 0;
	else
		return 1;
}

int check_satisfy2 (int N, int L, int pathway[L][N], int i, int mode, int mode1, int mode2, int irflag, int bdflag, int dgflag, int svflag, int pnflag, struct designability *designability, int Link_dist[N + 1], int degree_dist[2 * N + 1], struct solution *solution)
{
	int j, k;
	int state[L - 1][2 * (N + mode1) + 1];

	for (j = 0; j < L - 1; j++)
	{
		state[j][0] = state[j][N + mode1] = 1;
		for (k = 0; k < N; k++)
		{
			state[j][k + mode1] = pathway[j][k];
			state[j][k + N + 2 * mode1] = pathway[j + 1][k];
		}

		if (j != L - 2)
			state[j][2 * (N + mode1)] = pathway[j + 2][i];
		else
		{
			if (mode == 1)
				state[j][2 * (N + mode1)] = pathway[1][i];
			else
				state[j][2 * (N + mode1)] = pathway[j + 1][i];
		}
	}

	return check_satisfy1 (N + mode1, L - 1, state, i + mode1, mode, mode1, mode2, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, solution);
}

void combine_dist (int N1, int N2, int dist[N2 + 1], double Dist[N1 * N2 + 1], double Dist2[N1 * N2 + 1], int i)
{
	int j, k;

	if (i == 0)
	{
		for (j = 0; j <= N1 * N2; j++)
		{
			Dist[j] = 0.0;
		}
		for (j = 0; j <= N2; j++)
		{
			Dist[j] = dist[j];
		}
	}
	else
	{
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist2[j] = 0.0;
			if (j <= N2)
			{
				for (k = 0; k <= j; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
			else
			{
				for (k = 0; k <= N2; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
		}
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist[j] = Dist2[j];
		}
	}
}

int calc_design (int N, int L, int pathway[L][N], int mode, int mode1, int mode2, int irflag, int bdflag, int dgflag, int svflag, int pnflag, struct designability *designability, double Links_dist[N * N + 1], double Degree_dist[2 * N * N + 1], struct solution *solution)
{
	int i, j, Link_dist[N + 1], degree_dist[2 * N + 1];
	double Links_dist2[N * N + 1], Degree_dist2[2 * N * N + 1];

	designability->NS = 1.0;
	designability->min_links = 0;

	for (i = 0; i < N; i++)
	{
		if (mode == 2 && i == 0)
			continue;

		designability->min_link[i] = N;
		if (bdflag)
		{
			for (j = 0; j <= N; j++)
			{
				Link_dist[j] = 0;
			}
		}
		if (dgflag)
		{
			for (j = 0; j <= 2 * N; j++)
			{
				degree_dist[j] = 0;
			}
		}

		if (svflag)
		{
			solution[i].curr = 0;
		}

		if (check_satisfy2 (N, L, pathway, i, mode, mode1, mode2, irflag, bdflag, dgflag, svflag, pnflag, designability, Link_dist, degree_dist, &solution[i]) == 0)
			return 0;

		designability->NS *= designability->NS2[i];
		designability->min_links += designability->min_link[i];
		if ((i == 1 && mode == 2) || (i == 0 && mode != 2) || designability->NS2[i] < designability->NS2_min)
			designability->NS2_min = designability->NS2[i];

		if (bdflag)
		{
			if (mode == 2)
				combine_dist (N, N, Link_dist, Links_dist, Links_dist2, i - 1);
			else
				combine_dist (N, N, Link_dist, Links_dist, Links_dist2, i);
		}

		if (dgflag)
		{
			if (mode == 2)
				combine_dist (N, 2 * N, degree_dist, Degree_dist, Degree_dist2, i - 1);
			else
				combine_dist (N, 2 * N, degree_dist, Degree_dist, Degree_dist2, i);
		}
	}

	return 1;
}
